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ABSTRACT 



We search for an infrared signature of the transiting extrasolar planet 
HD 209458b during secondary echpse. Our method, which we call 'occultation 
spectroscopy,' searches for the disappearance and reappearance of weak spectral 
features due to the exoplanet as it passes behind the star and later reappears. 
We argue that at the longest infrared wavelengths, this technique becomes prefer- 
able to conventional 'transit spectroscopy'. We observed the system in the wing 
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of the strong i/3 band of methane near 3.6 /xm during two secondary echpses, 
using the VLT/ISAAC spectrometer at a spectral resolution of 3300. Our anal- 
ysis, which utilizes a model template spectrum, achieves sufficient precision to 
expect detection of the spectral structure predicted by an irradiated, low-opacity 
(cloudless), low- albedo, thcrmochemical equilibrium model for the exoplanct at- 
mosphere. However, our observations show no evidence for the presence of this 
spectrum from the exoplanet, with the statistical significance of the non-detection 
depending on the timing of the secondary eclipse, which depends on the assumed 
value for the orbital eccentricity. Our results reject certain specific models of the 
atmosphere of HD 209458b as inconsistent with our observations at the 3a level, 
given assumptions about the stellar and planetary parameters. 



1. INTRODUCTION 

The discovery of the transiting extrasolar planet HD 209458b (Charbonneau et al. 2000; 
Henry et al. 2000) has provided a unique opportunity to deduce the physical characteristics 
of one example of the so-called 'hot Jupiter' class of exoplanets. The transit geometry 
has allowed an accurate derivation of the stellar and planetary mass and radius (Brown 
et al. 2001b) and also provides an opportunity to measure the composition of the exoplanet 
atmosphere. The scale height of an atmosphere in hydrostatic equilibrium is proportional to 
T/g, where T is the atmospheric temperature and g is the surface gravity. Because stellar 
irradiation of the planet is intense at 0.046 AU (Seager & Sasselov 1998), the atmosphere 
could be as hot as T ~ 1400 K or more. Coupled with the relatively low surface gravity 
{g — 848 cm s~^), the atmospheric scale height is potentially as large as if ~ 750 km. 
Moreover, when viewed tangent to the limb, the atmosphere can be opaque over several 
scale heights at the wavelength of a strong absorption line. The effective absorbing area of 
the planet can thus be greater in a strong line by a detectable amount, producing a slightly 
deeper eclipse at wavelengths close to the line core (Seager & Sasselov 2000; Brown 2001; 
Hubbard et al. 2001). That is, the opacity of the exoplanct atmosphere will impose weak, 
but potentially detectable, absorption lines on the stellar spectrum as the planet passes in 
front of the star. Several investigators have attempted to detect this effect (Bundy & Marcy 
2000; Moutou et al. 2001). Using the Hubble Space Telescope to observe the primary echpse, 
Charbonneau et al. (2002) succeeded in detecting the sodium doublet in the atmosphere of 
HD 209458b. 



Based on observations collected at the European Southern Observatory, Chile (proposal no. 67.C-0196) 
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In addition to sodium, absorption in molecular features of water, methane, and carbon 
monoxide may also be observable (Hubbard et al. 2001; Brown 2001). Brown et al. (2002) 
attempted to detect absorption by carbon monoxide using ground-based observations. From 
their attempt and our own experience, we recognize that the terrestrial atmosphere can 
potentially frustrate such efforts, but we believe that good observing conditions and appro- 
priate observational and analytical techniques will enable successful measurements from the 
ground. The known timing of the transit (Brown et al. 2001b) is a boon to this technique, 
because it discriminates against false signals by providing a temporal modulation on the 
signal. 

In this paper we introduce a new observational approach involving the secondary eclipse, 
when the planet passes behind the star. The importance of the secondary echpse was first 
pointed out by Charbonneau et al. (2000), but no measurement attempt has been reported 
in the literature. We have performed moderate resolution 'occultation spectroscopy' to 
measure the modulation (i.e., disappearance and reappearance) of the potential exoplanet 
line spectrum as the planet is hidden and revealed by the star. We have also performed 
lower resolution spectroscopy during and near secondary eclipse using the NASA Infrared 
Telescope Facility (IRTF), attempting to isolate the broad infrared (IR) flux peaks from the 
planet's presumably hot atmosphere at wavelengths where the opacity is minimal. This paper 
introduces the concept of occultation spectroscopy for extrasolar planets and reports the Very 
Large Telescope (VLT) spectroscopic results; lower resolution occultation spectroscopy of the 
IR flux peaks using the IRTF observations will be the subject of a second paper. 



2. OCCULTATION SPECTROSCOPY 

In this section, wc briefly explain our rationale for attempting occultation spectroscopy 
at secondary eclipse, as opposed to more conventional transmission spectroscopy during 
primary eclipse. At visible and near-IR wavelengths, the self-emitted flux of the planet is 
negligible, making occultation spectroscopy of thermal emission impossible. However, many 
of the strongest molecular features of interest in exoplanets (e.g., the very strong fundamental 
1^3 band of methane) occur at longer IR wavelengths. With increasing wavelength, the total 
signal detectable using occultation spectroscopy becomes as large, or larger, than that from 
transit spectroscopy. To see this, consider a spectral line having sufficient strength to be 
opaque N scale heights higher in the atmosphere compared to the nearby continuum when 
viewed tangent to the limb. Cases of interest are for large iV, so we can safely assume that 
the line is also optically thick when viewed at normal incidence to the exoplanet atmosphere. 
Transit spectroscopy looks for this absorption against a stellar continuum fiux proportional 
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~ 27r RNHBl, where H is the atmospheric scale height of the exoplanet, 2'kRNH is the area 
of the annulus of the planetary atmosphere, and B\ is the Planck function for the stellar 
continuum temperature at wavelength A. Occultation spectroscopy, on the other hand, looks 
for the modulation of the exoplanet spectrum itself, with continuum flux proportional to 
~ TiR^B^, where is the Planck function for the exoplanet continuum temperature. With 
increasing wavelength, the increase in B^ and the large area of the exoplanet (vri?^) dominates 
the area of the atmospheric annulus {2t^RNH) in spite of the greater stellar intensity {B*^. 
Evaluating the ratio of these fluxes for A?" = 3, if = 750 km, R— 1.0 x 10^ km, and for stellar 
and exoplanet temperatures of 6000 K and 1200 K, respectively, we find that the ratio equals 
unity for A = 2.6 //m. A rigorous treatment would of course include the line profiles and 
many other factors, but the inevitable tendency is that, for wavelengths longward of the K 
photometric band, occultation spectroscopy becomes increasingly favorable over transmission 
spectroscopy in terms of the amplitude of the potential signal. 

Occultation spectroscopy is essentially direct IR spectroscopy of the exoplanet, and 
in principle this could be done at most orbital phases, and also done for planets which 
do not transit. However, it requires detection of the small-amplitude exoplanet spectrum 
(the precise nature of which is unknown), superposed on a strong stellar background. The 
measurement becomes much easier if the signal is modulated predictably. The disappearance 
of the exoplanet spectrum during secondary eclipse is therefore an essential part of our 
technique. Note, however, that modulation by variations in Doppler velocity could also be 
used for all hot Jupiter systems, even ones without transits (Wiedemann et al. 2001). Note 
also that the concept of occultation spectroscopy is not new; it was successfully used to 
obtain the IR spectra of Pluto and its moon Charon (Sawyer et al. 1987). 



3. OBSERVATIONS 

We observed two secondary eclipses, on UT 2001 July 8 and UT 2001 July 15, with the 
ISAAC spectrometer on the VLT (Antu) at Cerro Paranal in Chile. The observations were 
performed in visitor mode, with real-time decisions regarding nodding frequency, integration 
time, etc. Both nights were clear, with 10 — 20% relative humidity, and seeing in the range 
from 0.5 — 0.7 arcsec. The zenith column of precipitable water vapor, as measured from 
telluric water lines in our spectra, was ~ 6 mm on July 8 and 1.6 mm on July 15. 

ISAAC is capable of imaging and spectroscopy in the wavelength region between 1 
and 5 //m (Moorwood 1997). We used spectroscopic mode with a resolving power of 3300, 
covering the 3.5 to 3.75 /im range (2667 — 2852 cm'^). Both eclipses occurred within an 
hour of transit across the local meridian, where the star reached a minimum airmass of 1.38. 
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The spectra maintained good quality to the largest airmass we observed (2.51). 

At approximately half- hour intervals, we observed a comparison star, HD 210483. The 
comparison star is within a degree of HD 209458, with nearly the same B and V magni- 
tudes and photometric colors (Table 1). We observed both stars using a conventional 'nod' 
technique, placing them alternately at 'a' and 'b' positions on the slit, in the order 'abba', 
with 60-second integrations at each slit position. On July 8 we obtained 36 'abba' sets of 
HD 209458, and 20 'abba' sets of HD 210483, while on July 15 we recorded 41 and 24 sets, 
respectively. About once per hour, we obtained spectra of a continuum lamp for fiat-fielding. 
We did not rely on the standard fiats; instead, we recorded fiat field exposures immediately 
after each block of HD 209458 observations without moving the telescope or instrument, in 
order to avoid small but non-negligible flexure effects. 



4. MODELED SPECTRUM 

We do not expect to achieve sufficient precision to detect the planetary spectrum in 
each individual spectral resolution element; instead, we look for the candidate signal by 
performing a least-squares fit of the observations to a model 'template'. We computed the 
template as the emergent flux from a model atmosphere for the exoplanet, using a series of 
speciflc intensity spectra at different values of (the cosine of the zenith angle) . These were 
computed from the formal integral of the radiative transfer equation, using an LTE source 
function; we neglected all scattering terms, since the long-wavelength region is of interest 
and we are assuming a cloud-free model. Flux was computed by a quadrature integration of 
the intensity spectra. The spectrum synthesis was done at very high spectral resolution, and 
convolved to the 0.8 cm~^ resolution of our ISAAC data using a Gaussian instrument profile. 
We adopted the parameters for the star and planet given by Cody & Sasselov (2002). 

Model atmospheres for the planet can vary widely, with one major difference being 
the treatment of clouds and aerosols. The maximum thermal contrast in IR spectra is 
generally obtained from clear atmospheres; one reason for this is the fact that scattering 
processing in clouds tend to wash out strong absorption features. Moreover, the presence 
of clouds high in the atmosphere effectively limits the detection of methane by allowing 
only the small fraction of the column above the cloud to be probed by the observer. We 
have therefore used the most optimistic case of the clear atmosphere as a reference point. 
Our fiducial model follows the calculation of Goukenleuque et al. (2000), with the surface 
gravity appropriate for HD 209458b. This model includes irradiation by the star, but uses 
only gaseous opacities; condensates and particulates are assumed to have settled out of 
the atmosphere forming a deep cloud layer. Figure 1 compares the pressure/temperature 
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structure from our fiducial model with a recent model from Sara Seager (2002, private 
communication). The temperature vs. optical depth structures of these models are similar, 
but the fiducial model requires a much higher pressure (due to lower opacity) to attain the 
same optical depth as the Seager model. (We also inspected a model for HD 209458b kindly 
sent to us by David Sudarsky (2002, private communication); it is slightly hotter than the 
Seager model at all depths.) 

We calculated the number densities of methane in the fiducial model, using the simple 
analytic formulae given by Burrows & Sharp (1999). The methane mixing ratio peaks at 
~ 7 X 10~^ at 25 mbar pressure, similar to the values for the 51 Peg model (Goukcnicuque 
et al. 2000). In computing the template, we used the current best-available experimental 
methane line parameters from Brown, Dulick, & Devi (2001a), now incorporated into HI- 
TRAN (Rothman et al. 1998), and the calculation followed Wiedemann et al. (2001) in other 
respects (partition functions, etc.). Since line parameters are critical to calculation of the 
template, we calculated spectra of hot methane (773 K) at ISAAC resolution for laboratory 
conditions, using both our adopted line data and also theoretical line parameters from the 
spherical top data system theory (Wenger & Champion 1998). We compared these calcu- 
lated spectra to laboratory measurements of hot methane at 773 K, measured by Steyert & 
Renter (2003, private communication) at the same resolution. Based on this comparison, we 
have confidence in the HITRAN methane line data used to compute our model spectrum. 
Since our spectral region also exhibits absorption due to water, we experimented using water 
lines in our template, also based on HITRAN. Because the HITRAN database is incomplete 
at high temperatures, we also experimented with using HITEMP, a HITRAN line database 
extension for higher temperatures, to obtain water line data. This results in abundant weak 
absorption lines that do not significantly change the template. We have therefore neglected 
water absorption in the analysis until more complete information is available for individual 
hot water lines in this spectral region. 

We expressed the exoplanet flux spectrum as a ratio to the stellar fiux. We computed 
the stellar fiux from a gray atmosphere at T = 5800 K, reducing the adopted 6000 K stellar 
temperature by 200 K to allow for the decrease in brightness temperature from the visible 
to 3.6 fira, as documented for the solar spectrum by Vernazza et al. (1976). 

Figure 2 shows the modeled template spectrum, with and without water lines from HI- 
TRAN included. In the results quoted and discussed below, we have used the 'methane only' 
version of the model template, for reasons given above. However, our results and associated 
errors do not change significantly if 'HITRAN water' or 'HITEMP water' is included in the 
model template. 

We computed methane spectra from models intermediate between our fiducial model and 
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the Seager model. If we increase the continuous opacity in the fiducial model, while preserv- 
ing the same temperature-optical depth relation, the methane bands weaken rapidly (again, 
see Figure 2). The greater continuous opacity pushes all line formation to lower pressure 
layers, giving lower column densities 'above the continuum.' The mixing ratio of methane 
decreases at low pressures (Burrows & Sharp 1999), further lowering the methane column 
density. Our fiducial model is a limiting case of an exceptionally clear atmosphere, and 
many other models for the exoplanet will not exhibit significant methane features, thereby 
escaping the test of our present analysis. 

5. DATA ANALYSIS 

Analysis of the HD 209458 observations took place in four stages: 1) extraction of 
spectra from the 2-D data frames, 2) removal of the telluric (and the few intrinsically stel- 
lar) absorption features from the spectra, yielding residuals which potentially contain the 
exoplanet spectrum, 3) estimation of the 'amplitude' of the exoplanet's spectrum in each 
residual spectrum, and 4) fitting of the aggregate results to a secondary eclipse curve. Since 
our present (and future) results depend critically on the fidelity of our data analysis, we 
describe this process in some detail, and we also describe checks made to insure the integrity 
of the analysis. 

5.1. Spectral Extraction 

The spectral frames were first cleaned of 'hot pixels' and temporary anomalies such as 
energetic particle events. Temporal sequences of frames at a single nod position were exam- 
ined to locate hot pixels based on their variability, and these pixels were given zero weight in 
subsequent analysis. Energetic particle events were identified and repaired separately using 
a median filter applied to the temporal sequence. The cleaned 'abba' set was combined by 
adding the 'a' frames and subtracting the 'b' frames to produce a difference image, which 
was then divided by a fiat-field frame. The fiat-field frame was computed as an average of 
the individual fiats during the night. We verified that temporal variations in the fiats were 
negligible for purposes of our analysis. 

At an excellent IR site such as Paranal, changes in the thermal background are ap- 
proximately linear with time during the several minutes required to record an 'abba' set. 
Consequently, most of the background emission is removed by construction of the difference 
image a-b-b-|-a. However, we found that a second order correction was necessary to remove 
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the background completely. We compute the residual background spectrum by summing the 
region between the a and b slit positions in the difference image. The relatively long slit 
used in ISAAC facilitates the precision of this correction by providing a 17-arcsec length of 
background. This residual background spectrum is subtracted from each row of the two- 
dimensional difference image prior to the extraction of the spectra. Note that our technique 
differs from usual IR practice, which would be to subtract the (negative) 'b' spectrum from 
the 'a' spectrum after their extraction from the difference image, which would automatically 
subtract the background to second order. However, this would not maintain separate 'a' 
and 'b' spectra, which we desired in order to verify that our results are independent of the 
position of the spectra on the detector array. 

We extracted the 'a' and 'b' spectra from the difference image using the optimal extrac- 
tion formalism of Horne (1986). The spatial profile needed in the optimal extraction was 
constructed by fitting polynomials in the intensity vs. wavelength distribution at each slit 
position, and then sampling these curves in the spatial direction at each wavelength. In the 
process of extracting the spectra, we noted that the approximate signal-to-noise ratio of each 
spectrum (~ 200) was close to the value expected based on the noise from Poisson statistics 
of photoelectrons from the background. 

At this stage of the analysis, the spectra were separated into sets corresponding to 
July 8 or July 15, HD 209458 or HD 210483, and 'a' or 'b' slit positions. Each of the 
8 sets was analyzed independently of the other sets. The individual spectra in each set 
were shifted slightly using spline interpolation, so that the telluric absorption features were 
all coincident. Frequencies in wavenumber for the telluric lines were derived from high- 
resolution solar Fourier transform spectra (Livingston & Wallace 1991), convolved down to 
the resolution of the ISAAC data. (This accounted for any blending of telluric hnes, which 
can change their effective wavenumbers.) Wavenumbers for each point in our ISAAC spectra 
were derived by spline interpolation using the telluric lines as standards; the error of this 
procedure was ~ 0.1 cm~^, much less than the 0.8 cm~^ spectral resolution. A single accurate 
wavenumber scale was thereby inferred for all the spectra in each set. Sample spectra from 
both nights are shown in Figure 3. 



5.2. Telluric Correction 

The intensities in each spectrum fluctuate due to variable losses at the ISAAC slit, as well 
as changes in the telluric line and continuum absorption. For each set of spectra, we made an 
intensity normalization, followed by a correction to remove the telluric absorption. A pseudo- 
continuum was used as an aid to normalization; it was determined as a spline fitted to a 
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set of points in each spectrum having the greatest intensity in pre-determined wavenumber 
intervals. We then used a wavelength-independent scale factor to normalize each spectrum 
so that the integral intensity under the pseudo-continuum was constant for all spectra. That 
is, we enforce the condition that the total intensity under each pseudo-continuum is equal 
to the total intensity under the average pseudo-continuum; the normalization factor /j is 
calculated from 

f^-^ (1) 

where c represents the average pseudo-continuum, q represents the fitted pseudo-continuum 
to spectrum i, and fi is the factor by which each spectrum i is multiplied. Thus, the 
normalized spectra are calculated by 

rii = fiSi (2) 
where rii is the normalized spectrum and Sj is the original spectrum. 

In the normalized spectra, at each wavelength element, we fit a least-squares line to 
the natural log of intensity with airmass, lnnj(x. A). The fitted slope bx was used to correct 
each intensity value to the minimum airmass in that set of spectra (this corrects for airmass- 
dependent telluric absorption). That is, the airmass-corrected spectrum n[ is obtained from 

= exp{]nni{x,X) -bx{x - Xmin)) (3) 

Each set of corrected spectra was then averaged. Residual spectra were calculated by sub- 
tracting the average spectrum and dividing the difference by a continuum fitted to the 
average spectrum. This produced residual spectra given by 

n ■ - n' 
n = 

which presumably contain the candidate exoplanet signal. 

It is important to note several properties of the residual spectra. Division by the pseudo- 
continuum of the average spectrum leaves the residuals expressed in units of the stellar 
continuum. However, since we divided by only one continuum curve per data set, temporal 
changes in intensity were not masked by fluctuating errors in continuum flts. Also, the 
subtraction of the average spectrum means that the sum of the residuals at each wavelength 
is identically zero. Thus, only variations in the spectra survive this process. In particular, 
the exoplanet spectrum will survive (since it varies due to the eclipse), but the intrinsically 
stellar lines (Farmer 1989) are subtracted. 



(4) 
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5.3. Higher-order Corrections 

Several additional corrections were needed before the residuals could be compared with 
the model spectrum. First, we noticed that a broad absorption feature was evident in some 
of the residual spectra, and it corresponded to a blend of water lines near 2720 cm~^. The 
depth of this feature relative to the stellar continuum was ~ 30% greater on the more humid 
night (July 8), and variations of the depth of the feature on either night were about ±2%. 
To correct these variations, we obtained a low-noise residual water spectrum by differencing 
the natural logarithms of the average spectrum for July 8 and July 15 (for both stars, and 
for both the 'a' and 'b' spectra, separately). This exploits the (fortuitous) fact that the 
precipitable water vapor column was appreciably different on the two nights, but telluric 
methane absorption was the same. (Unlike the concentration of water in the troposphere, 
which is highly variable, methane is effectively constant because it is long-lived and therefore 
well-mixed.) We tested each residual spectrum to see whether a significant correlation existed 
with the water vapor spectrum. If it did, then we removed it from the residuals by least- 
squares fitting. 

We also noticed that the residuals sometimes exhibit variations at low spatial frequencies 
(see upper panel of Figure 4). These basehne effects (typically ~ 0.2 — 0.3%) are not 
surprising, because our analysis treats each wavelength point independently. (We do not 
introduce artificial baselines - e.g., by continuum fitting - but neither do we suppress real 
baseline variations.) A Fourier analysis also revealed a spike of excess power at the Nyquist 
frequency (0.5 cycles per pixel), which is commonly seen in IR data, usually indicating a 
difference in the two output amplifiers reading alternate columns of the detector array. We 
removed these effects by computing the Fourier transform for each set of residuals; then we 
zeroed-out an appropriate number of low frequency elements, as well as the DC offset (zero 
frequency) and the Nyquist frequency components. The number of low frequencies removed 
(up to 0.0079 cycles per pixel) was taken to be the maximum possible without impinging on 
frequencies where the exoplanet model spectrum exhibited significant power. The same filter 
operations were performed on the model template spectrum as on the observed residuals. 
Then we inverse-transformed the filtered residuals and template, which were then suitable 
for comparison. The effect of this Fourier filtering on the observed residuals is seen the lower 
panel of Figure 4. 

Even after the removal of variable water vapor absorption, and Fourier filtering, the 
residuals needed an additional correction. We discovered that they contained an excess of 
'outliers', i.e., points more than several standard deviations from zero. We exploited our 
comparison star observations to find those wavelengths which tended to become outliers 
and reject them from the analysis. The procedure was applied to the 'a' and 'b' spectra 
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separately, since we suspected low- amplitude 'hot pixels' as a likely cause of outliers. We 
applied a wavelength-dependent 'weighting mask' to the residuals, where the weight of the 
mask (normally unity) was set to zero when a particular wavelength was found to deviate 
from zero by more than 2.5(T in two or more spectra for both HD 209458 and HD 210483. 
The requirement that points are zero- weighted only if that wavelength tends to be discrepant 
in spectra from both stars assures that this procedure does not bias the results. About 7% 
of the points were zero-weighted using this method; we found that our final results were 
remarkably insensitive to the limits for the rejection. The reason for this is that the total 
number of outliers was much less than 7%, because a given point was always zero- weighted 
when it met the above criterion, even if it was 'well-behaved' in the majority of spectra. 
The standard deviation of the Fourier-filtered and masked residual spectra was in the range 
0.003 to 0.006. 



5.4. Fit to the Exoplanet Spectrum 

Prior to fitting the model template and Fourier filtering, the residuals are corrected for 
the heliocentric Doppler shift due to the relative motion of HD 209458 with respect to the 
Sun (Nidever et al. 2002), as well as the geocentric Doppler shift due to the relative motion 
of the Earth around the Sun. We have also corrected for the non-negligible Doppler shift due 
to the orbital motion of the exoplanet with respect to HD 209458. The residuals are shifted 
in wavenumber to the rest frame of the exoplanet, and the model spectrum (already in the 
exoplanet rest frame) is interpolated onto the same wavenumber grid. Note that although 
the correction is small, we have corrected each residual (each frame) for the Doppler effect, 
based on its time of observation. The residuals calculated for HD 210483 are analyzed in an 
identical manner. Although HD 210483 is not known to host a planet, we apply the same 
Doppler correction as we did to HD 209458, but we calculate the correction based on the 
time of observation of each HD 210483 residual. 

The comparison between each residual spectrum and the model exoplanet spectrum 
is made using linear least-squares. For this purpose we use the 'methane only' version of 
the model (see Section 4), since the state of the spectroscopic line data gives us the most 
confidence in this version. (Our overall results do not change with the 'methane plus water' 
version.) The result of a least-squares regression of the residuals versus the model spectrum 
gives a quantitative estimate of the degree to which each residual spectrum 'contains' the 
exoplanet model spectrum. We refer to these estimates as 'model amplitudes', since they 
represent the degree to which the model spectrum appears in each residual spectrum. The 
least-squares fit also estimates the random error in the model amplitudes, typically 1.5 — 2.5. 
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In other words, the signal-to-noise ratio for detecting the modeled exoplanet spectrum in 
a single set of residuals is about 0.5. This is sufficient precision to expect detection of the 
exoplanet spectrum in the average of ~ 50 residual spectra, provided that the signal-to-noise 
ratio increases as the square root of the number of spectra averaged. 

5.5. Checks on the Analysis 

We have checked our analysis procedures in several ways. Since we use linear least- 
squares, the analysis should have the property that averaging the best-fit model amplitudes 
for a set of residual spectra should give the same result as fitting to the average of those 
residual spectra. We verified that our results satisfy this identity exactly, if the time- variable 
Doppler shift applied to the residuals is neglected (with Doppler shifts included, the equality 
is approximate.) We have also confirmed that our numerical procedures do not attenuate 
potential signals. We added a synthetic signal, identical to the modeled spectrum, at the 
earliest practical stage of the analysis (immediately after the extraction of the spectra from 
the 2-D data frames). The least-squares solutions recover this signal (as a difference, out of 
eclipse minus in eclipse) with a best-fit model amplitude of near unity. 

Another important check on our analysis is to examine the nature of the noise, specifi- 
cally the distribution of noise as compared to a normal error distribution. This is important 
because our error estimates implicitly assume that the data reflect a normal error distribu- 
tion. Figure 5 shows the error distribution for all of the points (after masking) in the 'a' 
residuals of HD 209458 from July 8. Apart from an excess of points within la of zero, the 
distribution closely approximates a normal error curve. Similar distributions are found for 
the other sets of data. 

Figure 6 shows the distribution of 'model amplitudes.' In this case we construct the dis- 
tribution from the least-squares fits to all residual spectra for both HD 209458 and HD 210483 
from both nights; this comprehensive inclusion is needed to provide sufficient points to expect 
a reasonable approximation to the normal error curve. Considering the number of residual 
spectra represented (242), the data are in good agreement with the normal error curve. As 
a further check on the errors, we computed the variation in fitted model amplitudes which 
result from fitting Gaussian random noise whose standard deviation varied from 0.003 to 
0.006, the same as the range of our filtered and masked residuals. The results of this sim- 
ulation were in close agreement with the error estimates derived from the least-squares fits 
to the actual residual data. 



5.6. Secondsiry Eclipse Timing 



The timing of secondary eclipse depends on the orbital eccentricity. If the eccentricity is 
identically zero (as might be expected based on tidal circularization arguments) the secondary 
eclipse will occur exactly mid-way between primary eclipses. If the eccentricity is non-zero, 
it can occur up to ~ 90 minutes earlier or later (depending also on a;). Doppler observations 
to date (G. W. Marcy 2002, private communication) give e = 0.011 ± 0.015 with oj — 156°, 
which would imply — —31 minutes, i.e., a slightly earlier secondary eclipse. We computed 
the time of secondary eclipse as a function of the eccentricity, using P = 3.52474 days 
(G. W. Marcy 2002, private communication), and the zero-point from primary eclipse at 
To = 2451659.93675 HJD (Brown et al. 2001b), and we corrected for light travel time. We 
checked our calculations against a primary eclipse ephemeris by David Charbonneau (2002, 
http://www.astro.caltech.edu/~dc/frames.html), obtaining essentially identical results, and 
we did a similar check for the time difference introduced by using a non-zero eccentricity 
(D. Charbonneau 2002, private communication). Given a calculated time for secondary 
eclipse, we fit a simple eclipse curve to our model amplitude data by linear least-squares. 
Our eclipse curve uses two levels (in- and out-of-echpse) connected by straight lines. The 
time from first to fourth contact was taken to be 184.3 minutes, and from second to third 
contact 132.2 minutes (Charbonneau et al. 2002). Note that, due to the subtraction of an 
average spectrum in our data analysis, only the amplitude of the eclipse curve is significant; 
the 'zero-point' in the fitted eclipse curve is not meaningful. 



6. RESULTS AND DISCUSSION 

Given the precision achieved by our analysis, we would expect to detect the exoplanet 
spectrum if it is represented by the irradiated, cloudless, low-albedo atmosphere, with ther- 
mochemical equilibrium abundances of methane (and water), and also provided that the 
timing of the secondary eclipse is known and the stellar and exoplanet parameters (radii, 
etc.) are exactly as adopted. We first conclude that this precision can indeed be achieved 
using ground-based observations. The recent results and discussion of Brown et al. (2002) 
imply that ground-based detection is feasible, since these authors came within a factor of ~ 3 
of the required precision using data from only one night in poor weather. Our experience, 
using a somewhat different technique, confirms that the requisite precision can be obtained 
if good observing conditions prevail. 

Our results for the secondary eclipse are illustrated in Figures 7 and 8, showing the 
derived 'model amphtudes' for each HD 209458 spectrum ('a' and 'b' for both nights) ver- 
sus time from the assumed mid-point of secondary eclipse. Recall that the error bars for 
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individual spectra were derived from the least-squares fits to the residuals after telluric cor- 
rections, Fourier filtering, and masking. We fit an eclipse curve to the aggregate results at 
each assumed eclipse time. (We fold the curve about the mid-point of secondary eclipse for 
the figures.) Also, there is no evidence that the results from the 'a' and 'b' spectra analyzed 
separately are significantly different. Figure 7 shows the result for e = 0.011, i.e., taking the 
slightly non-zero eccentricity from the Doppler data at face value, again for all spectra. The 
solid hne is the least-squares estimate of the echpse amplitude, —0.1 ± 0.3, and the dashed 
line shows the levels corresponding to an eclipse amplitude of unity (i.e., if the model tem- 
plate spectrum correctly represents the exoplanet, and the timing is correct.) Figure 8 shows 
the corresponding result (0.5 ± 0.4) under the assumption that the eccentricity is identically 
zero. In each case the reduced chi-squared of the fit is ~ 1.5, indicating that the scatter of 
the data is only slightly larger than the independently-assigned error bars. Note that each 
plot for HD 209458 is accompanied by a similar fit to the comparison star, shown in the 
lower panels. In both cases, the comparison star shows no significant change, as expected. 

Adopting e = 0.011 from the Doppler data, our results for non-detection of the secondary 
eclipse (—0.1 ±0.3) exclude the model template spectrum at the > 3(T level. Given our error 
distributions, this case is firmly rejected. Assuming an eccentricity of zero, non-detection 
of the secondary eclipse is ambiguous (0.5 ± 0.4). However, even if we explore eccentricities 
over the entire plausible range from 0.0 to 0.03, and assume a; values in the first and second 
quadrants, the fitted secondary eclipse amplitudes never reach unity. Instead, the amplitudes 
vary from —0.2 to -1-0.9, and the amplitudes for the comparison star similarly range from 
—0.5 to -1-0.9. At a given echpse time, the amplitudes for the two stars are uncorrelated, as 
expected when sampling the noise envelopes resulting from fitting independent random data 
to eclipse curves of variable timings. We therefore conclude that the secondary eclipse is not 
detected using our 'methane template.' Given the similarity between HD 209458 and the 
comparison star in terms of the noise envelope, and noting the nearly identical reduced chi- 
squared values, we conclude that HD 209458b probably does not exhibit significant methane 
absorption features, in agreement with the models of Sudarsky et al. (2002). Methane 
absorption in the combined light spectrum (star plus planet) has a total equivalent width of 
no more than 0.025 cm~^ within our absorption bandpass (2667 — 2852 cm~^). 

Deviation of the exoplanet atmosphere from our fiducial model is not difficult to ex- 
plain. As discussed in Section 4, models of this exoplanet must be exceptionally clear if 
significant methane absorption is to appear in their spectrum. In this respect our fiducial 
model represents a limiting case. Our results are certainly consistent with the suggestion 
of a cloudy atmosphere based on the observed low sodium abundance (Charbonneau et al. 
2002). However, even in a clear atmosphere there may be other viable explanations for both 
the sodium result (Barman et al. 2002), and our present result (e.g., photochemical deple- 
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tion of methane). Nevertheless, it is significant that we can now begin to hmit the range of 
parameters of viable exoplanet models using ground-based observations. 
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Table 1. Stellar Photometric Data. 



Parameter 


HD 209458 


HD 210483 




7.648 


7.586 




0.594 


0.585 


b-y2 


0.361 


0.391 


ml2 


0.174 


0.175 


cl2 


0.362 


0.354 



^Hipparcos Catalogue (Ferryman et al. 
1997) 

^uvbyP Photometric Catalogue (Hauck & 
Mermilliod 1998) 
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Fig. 1. — Temperature vs. pressure structure for our fiducial model for HD 209458b (right- 
most curve), and for a model by Seager (2002). The asterisks mark the points of Rosseland 
optical depth unity (Seager model) and 3.6 /im optical depth unity (fiducial model). 
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Fig. 2. — Model template spectra containing only methane lines (lower solid curve), and both 
methane and water lines (dotted curve). The upper solid curve represents the 'methane only' 
template with a factor of three increase in the continuous opacity, but retaining the same 
temperature-optical depth relation. 
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Fig. 3. — Sample spectra of HD 209458 for both observing nights. Most of the hne structure 
is telluric, due to methane and water vapor, as noted by the line identifications. A few stellar 
hnes are also marked (Farmer 1989). Note that the telluric water vapor lines were stronger 
on 2001 July 8. The spectrum from 2001 July 15 has been shifted by +0.5 for plotting 
purposes. 
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Fig. 4. — Residual spectra corresponding to the example spectra shown in Figure 3. Upper 
panel shows the residuals from the extracted spectra after telluric correction. Lower panel 
shows the same residuals after the Fourier filtering process. The residual spectrum from 
July 15 (dashed line) has been shifted by -0.03 for plotting purposes. 
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Fig. 5. — Noise distribution of the residuals for the target 'a' spectra on the UT 2001 July 8. 
Solid hne indicates the theoretical value based on the Gaussian probability distribution. 
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Fig. 6. — Noise distribution for computed values of model amplitude. All measurements 
(242) are included: both stars, both nights, and for both the 'a' and 'b' spectra. Solid hne 
indicates the theoretical value based on the Gaussian probability distribution. 
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Fig. 7. — Fitted secondary eclipse curves (solid line) for HD 209458 (upper panel) and the 
comparison star (lower panel). These fits assume the current best-fit Doppler value of the 
orbital eccentricity, e = 0.01, resulting in the eclipse being 31.1 minutes earlier than that of 
an orbit with zero eccentricity. The dashed lines represent an eclipse amplitude of unity. 
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Fig. 8. — Fitted secondary eclipse curves (solid line) for HD 209458 (upper panel) and the 
comparison star (lower panel). These fits assume that the orbital eccentricity is zero. The 
dashed lines represent an eclipse amplitude of unity. 



